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Abstract 

We study nano-pattern formation in a stochastic model for adsorption-desorption processes with interacting 
adsorbate and hyperbolic transport caused by memory effects. It is shown that at early stages the system 
manifests pattern selection processes. Stationary stable patterns of nano-size are analyzed. It was found that 
multiplicative noise satisfying fluctuation-dissipation relation can induce re-entrant pattern formation related to 
non-equilibrium transitions. According to obtained Fokker-Planck equation kinetics of island sizes in a quasi- 
stationary limit is discussed. Analytical results are compared with computer simulations. 



1 Introduction 

From theoretical and experimental observations it is known that reaction-diffusion systems play an important role 
in the study of generic spatiotemporal behavior of non-equilibrium systems. Usually such models admit main 
contributions related to both local dynamics (chemical reactions type of birth-and-death processes) and mass 
transport. Novel experimental methods, such as field ion microscopy, scanning tunneling microscopy can be used 
to monitor chemical reactions on the metal surfaces with atomic resolution. 

In adsorption-desorption processes when material can be deposited from the gaseous phase such experimental 
methods allow one to investigate formation of clusters or islands of adsorbed molecules/atoms \T\. Such islands can 
have linear size of nanometer range [2]. In Refs. [SJjT] it was experimentally shown that nano-patterns on solid surface 
and nano-islands in adsorbed mono-atomic layers can be organized. It was found that lateral interactions between 
adsorbate play an important role in processes of pattern formation on metal surfaces. The adsorbate presence 
can modify the local crystallographic structures of the substrate's surface layer producing long range interactions 
between adsorbed atoms and their clusters (see for example, Refs. [8}|l0]). It was observed experimentally that 
nanometer-sized vacancy islands can be organized in a perfect triangular lattice when a single monolayer of Ag was 
exposed on Ru(OOOl) surface at room temperature [?]. The origin of such self-organization of vacancy islands is an 
elastic interaction between edges of neighboring islands [?]. Such systems can be widely used in electronic devices 



in industrial applications. Despite theoretical studies of nano-patterns were widely reported (see Refs. 11-16 and 
citations therein) no detailed descriptions of elastic effects leading to formation vacancy/adsorbate islands lattice 
were presented. 

It is known that such kind of processes can be studied by multiscale modeling including: ab-initio calculations 
allowing one to set parameters of primitive cell at zero temperature and some energetic parameters; molecular 
dynamics simulations monitoring microscopic processes of interactions of atoms; Monte-Carlo approach studying 
physics on bigger time scales. Unfortunately the corresponding computational time remains excessive and presented 
methods can be used to describe the system only on micron scales. To study pattern-forming transitions one can 
use mesoscopic approach to modelize self-organization processes on large sizes. It can be generalized by taking 
into account corresponding fluctuating sources representing an influence of a microscopic level or fluctuations of 
environment 
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Properties of pattern formation in systems of adsorption-desorption type with dissipative dynamics were studied 



previously 11 12 . Pattern formation in stochastic systems with internal multiplicative noise obeying fluctuation- 
dissipation relation was discussed in Refs. 13 16 , 17]. It was shown that internal multiplicative noise can sustain 
stationary patterns of nano-size range; the corresponding patterning processes are well related to non-equilibrium 



transitions induced by internal multiplicative noise (entropy driven phase transitions 18-22]) 



Recently modifications of this mesoscopic approach were proposed, where in the framework of continuous model 
one can describe processes at microscopic and diffusion scales. This approach is known in literature as modified 



phase field crystals method 23 -26 . It can be used to study crystal growth, solidification processes, motion of defects, 



phase transitions in crystalline systems (having special kind of order). Its main advantage is in introduction of a 
crystalline order with elastic interactions leading to formation of crystal lattice and of inertia effects related to a 
finite propagation speed 27-29 . Inertial contribution is essential at early stages in the system dynamics. On the 



diffusion time scale one can observe patterning with a formation of crystalline lattice regarding to the Hook's law. 
Formally, the elastic effects are modelized by expansion of a pair interatomic potential in the Fourier space up to the 



fourth order in wave-numbers [23 24 30 . As a result the spatial interactions are described by the Swift- Hohenberg 
operator {q^ + V^)^, V = d/dr, qq is a constant related to period of spatial patterns. Infiuence of both internal and 



external stochastic sources onto pattern-forming transitions in crystalline systems were discussed in Refs. 31 32 



It was shown that in such kind of systems due to finite propagation speed pattern selection processes are possible; 
the external noise can principally change topology of patterns realized in a pure deterministic case. 

For class of systems with adsorption-desorption and diffusion processes the typical deterministic equation, which 
will interest us is of the form 

dtx = fix) - V • J, (1) 

where x — x{r,t) is the local coverage at surface defined as the quotient between the number of adsorbed 
molecules/ atoms in a cell of the surface and the fixed number of available sites in each cell, a; < 1. The term 
f{x) stands for local dynamics and describes birth-and-death or adsorption-desorption processes; the flux J repre- 
sents the mass transport. 

Most of the theoretical studies deal with the standard Fick law J = —DWx, D is the diffusion constant. It 
is known that at / = the ordinary diffusion equation dtX = V • DVx has the unrealistic feature of infinitely 
fast (infinite) propagation. To avoid such unphysical situation the diffusion flux can be generalized by taking into 
account memory (non-equilibrium) effects 



M{t,t')D\/x{r,t')dt' 



(2) 



described by the memory kernel M{t, t') — t exp(— 1< — t'\/Tj). In such a case in the limit f{x) = one gets the 
flnite propagation speed ^DjTj. At tj —t^ the assymptotic M{t,f) = 5{t — t') leads to the Fick law J = —DVx 
with infinite propagation. As far as real systems (molecules, atoms) have flnite propagation speed one should use 
Eq.([2]) or the equivalent Maxwell-Cattaneo equation: Tjdt3 = —J — DVx |28 . Equations (lj2) can be combined 
into one equation of the form 



Tjdlx + Tjf'{x))dtx = fix) + V • Dix)Vx, 



(3) 



where prime denotes derivative with respect to the argument. As was pointed out in Ref. [28' this equation has 
some restrictions: (i) it typically does not preserve positivity of the solution a;(r,t); (ii) the damping coefflcient 
must be positive, i.e., fix) < tJ^ . 

In this study we are aimed to describe pattern formation processes in system given by Eq.(]3]) generalized 
by introducing multiplicative noise satisfying fluctuation-dissipation relation. We use suppositions of the phase 
field theory to construct effective model describing formation of nano-patterns with crystalline order and pattern 
selection processes. It will be shown that formation of stationary patterns in stochastic systems is well related to 
noise-induced effects, namely, non-equilibrium transitions. We shall argue that noise correlated in time can induce 
re-entrant ordering processes, where the stable or unstable patterns are realized in a window of control parameter 
values. Studying behavior of islands size as clusters of dense phase and vacancy islands of diluted phase we obtain 
its time asymptotics. Analytical results are verified by computer simulations. 

The paper is organized as follows. In Section 2 we describe the model of our system. Section 3 is devoted to 
study stability of patterns at early stages. In section 4 we derive the effective Fokker-Planck equation to study noise 
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induced effects in Section 5 and formation of stationary patterns in Section 6. In Section 7 we discuss evolution of 
islands size. We prove our analytical results by computer simulation in Section 8. Finally, conclusions are presented 
in last section. 



2 The model 



one assumes 



We consider a model where only one class of particles is possible. Following Refs. 11 13 16 33 -35 
that the particles can be adsrobed, desorbed, can diffuse and interact among themselves. Therefore, we introduce 
the scalar field describing dynamics of the local coverage at surface x{r,t) G [0, 1]. The reaction term incorporates 
adsorption and desorption terms as follows f{x) = kap{l — x) — kdX e:xjp{U{r)/T). Here fc^ and kd are adsorption 
and desorption rates, respectively; p is the partial pressure of the gaseous phase; U{r) is the interaction potential. 

The total flux J is a sum of both ordinary diffusion flux {—DqS/x) and flow of adsorbate (— (r'o/T)a;(l — x)'VU). 
Here the multiplier x{l — x) denotes that the flux is only possible to the (1 — x) free sites. Using the phase field 



crystals methodology 23 24 30] for the potential U one assumes 

U{r)/T = - J dr'u(r - r')a;(r') ~ -ex ~ fiV^x - ^cV^x, (4) 

where e = J dru{r), /i = (1/2T) J drr^u{r) and k = (1/4T) J drr'^u{r). To estimate e, fi and x one can 
suppose that there is maximal intensity u(0) = Umax and characteristic radius tq 12 , therefore, one has e ~ Umaxf^, 
fj, ~ erg, >c ~ /i?'^ The obtained expression can be rewritten in a more simplest form U{r)/T ~ —(l/4:e)CsHX, 
where Csh = (9o + V^)^ is the Swift-Hohenberg spatial coupling operator where Qq = e/x; here we take fJ. — 1 that 
results to relations: >c = l/4£, Qq = 4e^. 

Next, after substitution of all above terms we get a system of two equations kind of ([T] [2]) with 

f{x)=a{l-x)~xexp{~ex), a = kap/kd, (5) 

where the total diffusion flux J satisfies the following relaxation equation 

Tjdt3 = -J-Do[Vx-G{x)VCsHx]; G(a;) = ^a;(l - x). (6) 

The equivalent equation for the coverage takes the form 

Tjd^tx + ^ix)dtx = fix) + DoV • [Va: - Gix)\/CsHx], (7) 

where j{x) = 1 — Tjf'{x). In our stochastic analysis we are based on this equation generalized by introduction 
of a multiplicative noise in ad hoc form using the fluctuation dissipation relation. Therefore, the corresponding 
Langevin dynamics is described by the equation 

Tjdlx + j{x)dtx = ^{x; V) + g{x)C{r, t). (8) 

Here (p{x; V) = f (x) + DqW ■ [Wx~G{x)WCshx]; "y{x) = g'^{x) in order to guarantee that the fluctuation-dissipation 
relation is satisfied [36| ; ^ represents Langevin force. From the mathematical viewpoint a generalized Langevin 
equation defined with a memory kernel has a noise fulfilling the fluctuation-dissipation relation with the noise time- 
correlator as the memory kernel has. Therefore, in our case the memory kernel is reduced to the delta-function and 
hence, the process C is a white noise in time (with a zero correlation time, Tc = 0). It is known that mathematically 
defined white noise (with delta-correlation function) has no physical meaning [37| , for example: every next value for 
^ does not depend on it previous value; a spread in values of the noise C and a total power of the noise are infinite. 
Therefore, one needs to consider physical fluctuation source having small but flnite correlation time ^ 1 (see, 
for example [38| ). Hence, in our study we consider a generalized model of the Gaussian noise with zero mean and 
a correlator 

(C(r,t)C(r',0) = — exp (J-^) S{r-r'), (9) 

where is the noise intensity. Hence, a set of main functions defining the system behavior is reduced to the 
reaction term f{x) and the field-dependent diffusion coefficient G{x). 



^Here we assume that ro Ldif> where the diffusion length is -L^i/ = \/~DoJkd. For ro one has estimation ro ~ Inm, whereas 
Liif ~ l^lm [33]|34]. 
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3 Linear stability analysis 



Let us study fluctuation effects at early stages. The stability analysis for the system can be done considering a 
behavior of the averaged quantity (x) in the vicinity of the state xq related to unique minimum position of a bare 
potential V{x) ~ ~ dx' f{x'). Averaging Eq.ljsj) one gets 

(10) 



Tjdftix) + {jix)dtx) = (fix) + DoV ■ [Vx - G{x)\/CsHx]) + {g{x)Cir,t)) 



Making use of the Novikov theorem 39 , for the noise correlator we have {g{x)C{j^, t)) = (cr^/2) {j'{x)). Therefore, 
in the vicinity of the state a;o for the Fourier component ((5a;k(w)) of the quantity (dx) — {x{r, t)) —xq the dispersion 
relation takes the form 



- _ i-f{xa)Lj = f'{xo) - Dok'^il - G{xo)[qi 
It leads to the following dependencies for dispersion curves: 



Uj{k)-_ 



17(^0) 

2tj 



DoP{l~GixoM^P]^)-f'{xo) 



Tj 



4rJ 



-,1/2 



(11) 



(12) 



One can see that uj{k) can have real and imaginary parts, i.e. uj{k) = ^u){k) ± &Lu(k). The component 9fia;(fc) is 
responsible for oscillatory solutions, whereas ^Ld{k) describes a stability of the solution ((5xk(w)). Analysis of both 
3?a;(fc) and Quj{k) allows us to set a threshold for a wave-number where oscillatory solutions are possible, from one 
hand, and, from another one, to define a wave-number for the first unstable solution. From the obtained dispersion 
relations it follows that at fc = fco satisfying equation 



kiil ~ GixoM - k'„Y) = 



1 



fixo)-—f M + 



(13) 



two branches of the dispersion relation degenerate. The unstable mode appears at k = k^ obtained from the 
equation 

klil-GixoM-kl]^)^-^ 



f"'M 



(14) 



From Eq.(ll) it follows that a sign of the effective control parameter e^f = /'(xq) — {rja /2)f"' {xq) defines 
stability of the system. Moreover, as far as the sign of e^f depends on the reaction term derivatives in the vicinity of 
the state xq, one can conclude that noise-induced effects in the linear regime are possible at tj 7^ and f"'{xo) ^ 0. 

To make a quantitative analysis we consider firstly stability of the deterministic system setting — in the 
vicinity of the unique minimum position xq of V{x). As far as go increases with the lateral interaction energy, 
next, we take small values for e and study two cases related to small and large values for the adsorption rate. 
The corresponding dependencies of ^Lo{k) and 5a;(fc) are shown in Figjl] From Figjl^ obtained in the vicinity 
of diluted phase characterized by — modes with wave-numbers k > ko manifest decaying oscillations (here 
^uj+{k) = — 5Ra;_(fc) 7^ 0, 5a;+(fc) < 0). Inside both intervals k G [0,fcci] and k G [fcc2,fco] the corresponding 
modes are stable without oscillations (here ^uj±{k) = 0, ^uj±{k) < 0). In the domain k £ {kci,kc2) unstable 
modes characterized by 3a;+(A:) > and ^uj±{k) = appear. In the vicinity of the dense phase with — > 1 
(see FigjljD) there are no oscillatory solutions due to diuj±{k) — 0, but unstable modes can be realized. Here at 
small wave-numbers one has '^Lu±{k < kc) < 0, whereas ai k > kc all modes are unstable in the linear regime with 
Qa;+(fc) > In the stochastic case (cr^ 7^ 0) in the vicinity of the diluted phase the noise action results in a 
decrease of the number of oscillating modes shifting fco toward tt; the interval for unstable modes [feci, fcc2] enlarges 
(see thick lines in Figjlji) . In the vicinity of the dense phase the noise shifts position of the critical wave-number 
kc toward small k increasing the number of unstable modes (see thick lines in FigjljD). 

In Figure [2] we plot stability diagram a{k)) at fixed e for both deterministic and stochastic cases. Here above 
solid lines (thin for the noiseless case and thick for the stochastic one) the corresponding modes with k £ [fed, kc2] 
are unstable. Below dashed lines the stable modes manifest oscillating behavior. The first unstable mode fcco 

^In the simplest case of x = in expansion (j4| one can find that in the vicinity of the diluted phase with above set of parameters 
e and a all modes manifest decaying oscillations, there are no unstable modes, here ^u}±{k) < 0, Scj+(A:) > 0. In the vicinity of the 
dense phase one has ^uj± < 0, whereas (fc) > at fc > feg 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.1 0.2 0.3 0.4 

a) k b) k 

Figure 1: Real and imaginary parts of the dispersion relation uj{k) for deterministic and stochastic systems (thin 
and thick lines, respectively) at e = 2.5, Dq = 1, tj — 0.3: a) a = 0.03323, xq = 0.035 (thin and thick lines relate 
to cr^ = and = 0.2, respectively); b) a = 0.4055, xq = 0.75 (thin and thick lines relate to = and = 0.1, 
respectively) 
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Figure 2: Dependencies a{k) at cr^ = (thin lines) and cr^ = 0.2 (thick lines). Other parameters are: e = 2.5, 
Do = 1, Tj = 0.3. Domains of unstable modes and stable oscillating modes are bounded by solid and dashed lines, 
respectively. 
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Figure 3: Dependencies km{o:) (see left ordinate axis) as dashed lines starting at fc = kco shown for cr^ = 0.2 and 
threshold values as dependencies e{a) (see right ordinate axis) at tr^ = and = 0.2 (thin and thick solid lines, 
respectively). Other parameters are: Dq — 1, tj = 0.3. 




is realized when equation 5a;(fc)-|_ = has unique nonzero solution fcco- At small lateral interaction strength of 
the adsorbate and small adsroption rate only stable modes without oscillations and stable oscillating modes can 
be formed. Above the threshold ath = a{kco) unstable modes appear leading to spatial self-organization of the 
adsorbate. When the noise intensity increases above pattern forming processes can be observed at lower value for 
the adsorption rate; the corresponding patterns are characterized by large wave-lengths then in the noiseless case. 
The thick dotted line denotes the dependence of the wave-number fc,„ related to a position of the maximum for 
9a;(fc)-|_. At fc > kcQ the corresponding value km defines the most unstable mode giving a period of the spatial 
structure. The most unstable mode at fc > fcco is characterized by 

(15) 

It follows that k„i depends only on the adsorption rate and the lateral interaction energy. The noise action leads 
to a change in the position of the point fcco lying in the line a(fcm)- 

Dependencies of the period of patterns at different values for lateral interaction energy are shown in Fig|3] as 
dashed lines fcm(a) starting at fc = fc^o for c'^ — 0.2. Here choosing a fixed value for e from dependencies e(a) one 
can find both the threshold value ath £^nd the corresponding fc^o. Above the obtained fcco the related period of 
patterns A = 2i: /km can be found at fixed system parameters. It follows that as the noise intensity increases the 
unstable modes appear at lower values for both a and e. At elevated e the selected patterns are characterized by 
small period. 

As far as the first statistical moment illustrates oscillatory behavior one should expect that processes of pattern 
selection can be observed in such system. Considering properties of pattern selection we need to find dynamical equa- 
tion for the structure function S{k.,t) as a Fourier transform of the two point correlation function {5x{r,t)5x{Y',t)) 
and to study its oscillatory behavior at small times. To that end we obtain the linearized evolution equation for 
the Fourier components 5x\^{t) and 5x-\^{t) and compute S'(k, t) = {5xi!^{t)5x-\s_{t)) . Using Novikov's theorem, and 
holding only leading terms in it, one gets a dynamical equation 



Tjdls{\^,t) + ^{xo)dtS{]^,t) 

= 2 (^/'(xo) - ^/'"(xo) - D^e{l - G{x,)[ql - ef )) ^(k, t) + 2^27(^0). 
Its analytical solution can be found in the form 5(k, t) — Sq (x exp(— iTi7(k)t), where 



(16) 



i7(a^o) , 2{D^k^{l~G{x^)[ql~k^f)- f'{xo)) , f{x^ 

w(fc)± = ^ ± W + o^rixo) - (17) 
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As in the previous case 9t37 is responsible for stability of the system, whereas SR-cu relates to pattern selection 



processes. Unstable modes are defined according to Eq.(f4). The quantity ko is determined as solution of the 
equation 



1 



/'M-^/"'(xo) + 



Sr., 



(18) 



An evolution of the structure function at early stages is shown in Fig|4^. It is seen that during the system 
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= 0.033 c = 2.5, n =0.02 

= 0.033 c = 2.4, = 0.02 

= 0.050 c = 2.5, = 0.02 

= 0.033 c = 2.5, a' = 0.20 




b) 

Figure 4: (a) The structure function dynamics at early stages t = 1, 2.4, 3 (the time is increased from top to 
bottom) in the vicinity of the state = 0.035 at £ = 2.5, a = 0.03323, Dq = 1, ct^ = 0.02. (b) The structure 
function a.t t — 2.4 and different values for e, a and a^. 



evolution the main peak related to wave-number of most unstable mode is formed; its height is increased. At small 
time interval except main unstable mode having largest peak in S{k,t) additional peaks a.t k > ko correspond to 
selecting modes f32||40 . As far as these modes are stable during the system evolution they disappear. Therefore, due 
to hyperbolic transport pattern selection processes can be realized in the adsorption-desorption systems. In Figjljj 
we plot the structure function at early stages at different values for lateral interaction energy and noise intensity. 
As diagram in Fig|3] shows an increase in both e and a shifts the main peak position of S{k) toward large values 
for the wave-number. The noise action does not shift the main peak position. The peak spreads due to fluctuation 
effects resulting in more diffuse interfaces between two phases. Moreover, as the lateral interaction energy increases 
the pattern selection processes occur faster than at small e. The same is observed for large adsorption rate and 
noise intensity. 



4 The effective Fokker-Planck equation and its stationary solution 

To consider statistical properties of the system at large time scales or in stationary case one can use the corresponding 
Fokker-Planck equation. To obtain it we exploit the standard technique and use the equivalent system of two 
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differential equations (Tj6) written for both generalized coordinate Xi and momentum pi on a grid of the mesh size 
£: 



TjXi = Pi, 

7(x) 



Pi + 



Tj 



Pi = ^{xi;Vij/£) + g{xi)Ci{t), 



(19) 



here i enumerates cells. Then, the corresponding probability density function is given by the average of the density 
functional p{{xi\,{pi\,t) over noise: P{{xi\,{pi\,t) = {p{{xi\,{pi\,t)), where means average over initial 
conditions. To obtain an equation for the macroscopic density functional P we exploit the conventional device and 
proceed from the continuity equation: 



dp 
dt 



E 



_d_ 
dxi 



d_ 
dpi 



p = 0. 



Inserting Xi and pi from Eq. (19), we obtain 



dp 

dt 



C + E p. 



(20) 



(21) 



where the operators £ — Ci and ^ — J2i^iCi ^^'^ defined as follows: 

Pi d ^ d f _ ^ 



TD dxi dpi 
d 

dpi 



(22) 



Within the interaction representation, the density functional reads p 



as 



d_ 
di 



p, that allows us to rewrite Eq. (21 ) 



(23) 



The well- kno wn cumulant expansion method serves as a standard and effective device to solve such a stochastic 
equation 41 . Neglecting terms of the order 0{7i?), we get the kinetic equation for the averaged quantity {p{t)) in 
the form 



d_ 
di 



{P){t) 



(7^(^)7^(^'))d^' 



(p)W- 



(24) 



/< 



Integrating the right hand side, one gets 

\n{t)1i{t'))At' = j dr^^MOe^^^^AA.Oe"''" 

= f dr V C{T)N,e^^N,e-^\ 
Jo 

Within the original representation, the equation for the probability density reads 

t 



(25) 



d 



P = IC + a^ / dT^C(T) 



p. 



(26) 



Because the physical time is much larger than the correlation scale, one can replace the upper limit of the integration 
with oo. 

To proceed we use the procedure proposed in Ref. 42 ■ 44 to obtain the effective Fokker-Planck equation for the 

(27) 



stochastic model. Then, expanding exponents, we obtain for the perturbation expansion 

dP 



dt 



= C + C]P 
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where collision operator C is defined as follows: 



-(0) 



(28) 



here C 



4n) 



{[A, B] is the commutator); moments of the noise temporal correlation function C(r) are: 



= ((72/n!)/Q°°T"C(r)dr. 

To perform the following calculations we shall restrict ourselves to considering overdamped systems where the 
variation scales ts, i, Xs, Vg, 7s, /s, Ds, and gs of the time t, the coordinate r, the quantity x, the velocity v = p/tj, 
the damping coefficient 7(2::), the force f{x) and the noise amplitude g{x), respectively, obey the following conditions: 



fsts 9sts ^ ^ 



D gX s 



= e-' » 1. 



(29) 



These conditions means a hierarchy of the damping and the deterministic/stochastic forces characterized by relations 

fs gs 



e « 1, = e « 1. 



As a result, the dimensionless system of equations (19) takes the form 

dxi _i 



Respectively, the Fokker-Planck equation (27) reads 



where the operator 
has the components 



i i 

The collision operator is defined through the expressions type of ([28]): 

oo 



(30) 
(31) 

(32) 

(33) 
(34) 



(35) 



dvi 



After suppressing the factor e 2, the collision operator written with accuracy up to the first order in e <C 1 takes 
on the explicit form up to second order in e: 



(0) 



92 



5,; 



92 1 fdg 



dvidvi 



d_ 

dvi 



92 



dvjdvi 



(36) 



dxidvi gi \dxi 

Note that in further consideration we are interested in the behavior of the particular distribution V{{xi}, t), not 
the total one, P{{xi},{vi},t). The reduced distribution can be obtained according to the moments introduced as 
follows 

" (37) 



7'(")({a;J,0^ J Pi{x,},{v^},t)l[Kdv^], 
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where the integration is provided over set {vi}. Then, perfor ming corresponding manipulations with Eq. (27), we 
arrive at the recursive relations for the moments 7^*^"^ ({a;^}, t) 43 : 

-g-pin+i) 



dt 



dxi 



nip 



•p("-i) 



dxi 



(„-l) 



(38) 



Here V = 7^*^°^ ({x^}, t), the first moment 'P'^^\{xi\,t) can be considered as a flux of the probability density V, i.e. 
V^^^ = J . Indeed, taking the zeroth moment, we obtain the expected continuity equation 



(39) 



The first moment calculation leads to 



(40) 



To evaluate the second moment T'^^^ we put n = 2 and take into account only terms of zeroth order in e: 



\ It 

Combining all obtained terms we arrive at the system of two equations of the form 



(41) 



(42) 



In the stationary case with no probability density current one gets the solution of the quasi-Gibbs form 

Ues[x] 



Vlx] (X exp 



where the effective energy functional is 



U,f[x] ^~Jdv(^J dx'^{x';W) + ^7(2:)) 



(43) 



(44) 



When we know the effective functional Uef [x] or its general construction in quadratures it can be used to describe 
possible phase transitions in the system or to analyze possibility of pattern formation. 

5 Noise- induced transitions 



Let us consider non-equilibrium (noise-induced) transitions in a homogeneous system, setting Dq = 0. The analysis 
of the effective potential Uef{x) 



Jdx'f{x') ~ ^^fl^7(x) resufl s in noise-induced scenario for modality change 



of the stationary distribution function P{x). The corresponding diagrams for noise-induced transitions are shown 
in FigjSj Dashed lines in FigjS^ correspond to coexistence of two phases: equivalence of two minima of the function 
V{x) — — J dx' f{x') (thin line) and Uef{x) (thick line). The dash-dotted lines starting at the point B relate to the 
system parameters where the minimum of Uef at a; = Xm equals Uef{x = bo), where &o = is the lower natural 
boundary of the diffusion process x{t). From the mathematical viewpoint at such choice of the system parameters 
one has a unique minimum at a; = x^, but due to existence of boundary 60 one has an "effective minimum" of U^f 
at X = 60. This minimum relates to the probability density of the adsorbate absence. Therefore, the dash-dotted 
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Figure 5: Phase and bifurcation diagrams for non-equilibrium transitions at tj = 0.3: a) phase diagram in the 
plane (a,e) at pure deterministic case (thin hnes) and A^*-^' = 0.1, 0.5 (thick lines); b) most probably stationary 
values Xm versus a at different e and A^^^) = 0.1. Domains / and /// correspond to diluted and dense states, 
respectively; in the domain // both states are possible. 



lines mean an equivalence of both the absence of the adsorbate and the existence of the dense state with Xm 7^ 0. 
The dotted line denotes values for e and a at which x^ takes zeroth values and belongs to the boundary 60 ■ At large 
noise intensity [M^^^ — 0.5) the dash-dot-dot line starting at the point C denotes equivalence of two "effective" 
minima of Uef at x = &o a-nd at a; = 61, where &i = 1 is the upper boundary of the diffusion process. Therefore, the 
dash-dot-dot line corresponds to equivalence of both states absence of the adsorbate and pure dense state. From 
Figj5^ it follows that the colored noise results in a decrease of the critical values for the system parameters a and 
e (cf. cases of white and colored noise influence when A^^-'^^ = and ^A'^^^ ^ 0, respectively). Therefore, if the 
system is subjected to the colored noise, then non-equilibrium unimodal/bimodal transitions are possible at lower 
lateral interaction energies e and adsorption rates a. Bifurcation diagram illustrating stationary values Xm versus 
the adsorption rate a is shown in FigjHjD. Here for the solid line quantities ai and a2 belong to spinodals (solid 
lines in Fig[5^); uq corresponds to the coexistence line; aob (see dashed line) relates to the coexistence line of two 
states Xm 7^ and Xm ^o- 



6 Analysis of stationary patterns 

The extreme of lA^f [x\ relates to the stationary patterns Xst which can be computed setting the first variation of 



Uef[x\ with respect to x equal to zero, i.e. 5Uef[x\ — 0. From Eq.(44| it follows that stationary patterns can be 
obtained as solutions of the problem 

f{x) + DoV • [Vx ~ G{x)\/CsHx] - -^-^ f"{x) = 0. (45) 

In order to study stability of possible solution in the vicinity of the fixed point x^ related to extrema positions 
of the homogeneous distribution we assume x{r) — x^ oc exp(i>!r7'), where >f can have real and imaginary parts 



responsible for spatial modulations and stability of the solution, respectively. After linearization of Eq.(45) one has 



Do>i'{l - G{x^)[ql ~ >^2]2) ^ ^,(^^^^) _ (46) 

where we put cp = 1, At'"^^ = TcCT^. 

According to obtained values for >r from Eq.(46l we have computed a stability diagram shown in Fig {§1 for 



both stochastic system with Ai'-^^ = 0.1 and deterministic one. It is seen that in the case of the noiseless system 
there is unique line of a hyperbolic form dividing domains of unstable (U) an stable (5) patterns (see thin line). 
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Figure 6: (a) Stability diagram for stationary patterns: U and S denote domains of unstable and stable patterns, (b) 
Profiles of stationary patterns related to points A{a = 0.14, e — 3.0), B{a = 0.2, e — 3.0) and C{a — 0.2, e = 5.0) 
in plot (a) at Dq = 3.0, M'^^'> = 0.1. 



At large e it tends to the axis a = 0. It follows that increasing the lateral interaction energy or the adsorption 
rate independently one can induce formation of stable stationary patterns. If colored fluctuations (A^^^-' ^ 0) are 
not negligibly small, then the picture of stationary pattern formation can be changed crucially. Here the thick 
curve dividing domains of stable and unstable stationary solutions at large e tends to the dependence e{a) denoting 
formation of the first nonzero Xm (line (OB) in Figjs^). It follows that as lateral interaction energy increases at 
fixed small values for adsorption rate one can move from unstable domain into domain of stable stationary patterns. 
At further increase in e we move again into domain of unstable solutions. Therefore, a some kind of re-entrant 
pattern-forming transition is realized. Starting from domain of unstable stationary patterns at elevated values for 
a we move into the domain of stable stationary solutions. Continuing an increase in e we get a re-entrantce, where 
unstable stationary patterns can be observed. The corresponding part of the stability diagram denoting re-entrant 
formation of unstable stationary patterns (a "beak" ) is a part of phase diagram describing noise-induced transitions 
(see Fig[5^). At further increase in e one can move again from domain of stable stationary patterns into domain of 
unstable stationary patterns. Finally, at large a the re-entrant ordering (formation of stable stationary patterns) 
can be observed. One dimensional stationary patterns are shown in Figj6]D where Lmax is the maximal period for 
one from three patterns obtained at a and e related to coordinates of points A, B and C in Fig|6^. It is seen that at 
small a one should expect a formation of islands of adsorbate. Such situation is possible for deposition of molecules 
(polymers) . At large a relevant for deposition of atoms on metal surface islands of diluted phase (vacancy islands) 
can be organized. 



7 Quasi-stationary limit 



Let us consider evolution of islands of adsorbate or vacancies (inclusions) having radius R taking into account that 
the most probable state Xm is reached and the system fluctuates around it. We assume {xm — Xco) < Xm where Xco 
plays a role of thermodynamics concentration value, in other words a supersaturation is small. 

Let us introduce a quantity y = Xm ~ xr measuring a deviation of the inclusion form from a sphere of the 
radius R. Using the standard approach j45| one can assume the concentration in the vicinity of a spherical inclusion 
as follows: xr = Xoo(l + Rs/R), Rg = 2i/f2/T is determined trough the surface tension v, atomic volume fl and 
temperature T. Then, the diffusion flux j{y) in the vicinity of the inclusion coincides with velocity of its radius 
change: dR/dt = j{y). Using the standard formalism of the coalescence theory 45 , we can replace the derivative 
operator V in the flux j{y) = Dq\7{1 — G{xm)^ CsH{<lo,'^))y by R~^ . As a result we obtain an evolution equation 
for the radius of inclusion in the form 



^.1(1 
dt' R 



G{x„ 



R 



1 

Rr 



(47) 
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where t' = tDoY., S = 2vQ,Xoo/T, A = Xm — Xoo > plays the role of a supersaturation, Rc = S/A > is the 
critical radius of the inclusion, CsH{q, R^^) = {q^ + R^^) ; next we drop prime for the time variable. Equation 



(47) can be rewritten in the form R — —dW/dR, where a potential W{R) can be used in a further study. 

At first let us consider the simplest case of A = const where the critical radius Rc remains constant. In the 



stationary case we put i? = in Eq.(47). Such assumption relates to autonomous regime of islands formations. As 
far as A ~ 10"^ the critical radius is large but finite. One of the stationary solution having physical meaning is 
Rq = L, where L — >■ oo is the thermodynamic limit of the system size. There is another solution Rq responsible 
for formation of stationary patterns characterized by 3?>f(ck,e, 7^ 0. It can be computed from the condition 

Rq = G{xm){^ + 'f'RoY- The corresponding dependencies i?o(a) at Rq < Rc at different s and M^^^ are shown 
in Figj7|i. Its values are bounded by stability criteria 9>f(a, e, Al*^^^) < 0. It follows that Rq relate to large stable 
values for x^j^ Using the potential W{R) one can find that W{R) has two asymptotics: W{R 0) — !• oo, 
W{R — oo) — — oo. Therefore, the quantity Rc corresponds to a maximum of W, whereas Rq defines its minimum 
position. Hence, at R > Rc all previously organized inclusions at early stages start to grow; contrary, at R < Rc the 
corresponding patterns relax to spherical inclusions having fixed radius Rq. In order to illustrate the corresponding 
time dependencies of the inclusion radius we solve numerically Eq.(47) (see FigjTjD) 



As far as Rq has unique dependence versus the system parameters, one can solve Eq.(47) analytically rewriting 
it for the quantity p = R/ Rc in the form p = [p — po){p ^ 1): where po = Rq/Rc, time derivative are taken 
with respect to t' = t/Rc (next we drop prime for the time variable). Dividing variables, from the formal solution 



of Eq. ( 47 1 following asymptotics can be found: in the vicinity ot p — pc = 1 there are exponentially increasing 
deviations from the critical radius; in the vicinity of po the exponential decaying p(t) is observed. 

Let us analyze the system behavior considering time dependence of the supersaturation A{t) > 0. Hence there 
is a time dependent variable Rc{t). Next it will be convenient to move to new variables: z{t) = Rc{t)/Rc{0), u(t) = 
R{t)/Rc{t). As a new time variable we introduce r = 41nz. It yields the derivative du'^/dr — —u* + '-f{T)R'^dR/dt, 



where 7(r) = {R'^{0)z^dz/dt) ^. Using Eq.(47l and taking into account that there is unique value Rq, one has 
R'^dR/dt — {R— Ro){R/Rc — 1). Therefore, an equation for the variable u takes the form 

du'* 

— = -u'^ - Roj{u-l)+jRcu{u-l). (48) 
dr 

Let us estimate time asymptotics for both R{t) and Rc{t). To this end we consider the limit t ^ 1. It allows 
one to put dw'*/dT = or = — i?o7(u — 1) + jRcu{u — 1). As far as the term in the left hand side has no time 
dependent coefficient and u(r — > oo) = const it means that all coefficients in the right hand side should be constants 
at r — >■ oo. When we assume 7(r — >■ oo) = const the coefficient jRc should take constant value too. Unfortunately, 
from the definition of 7(t) we arrive that Rc{t) is an increasing function in time leading to non-physical situation 



^In the simplest case of Xm = we arrive at classical results of the coalescence theory [45]. 
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Figure 8: Snapshots of the system evolution (a) and evolution of the spherically averaged structure function (b) at 
a = 0.16. Insertion in plot (b) shows structure function at t = 100 and different values for adsorption rate. Other 
parameters are: Dq = 3.0, tj = 0.3, A^^^^ = 0.1. 



where all coefficients from the right hand side of Eq.(48) are constants, not time-dependent functions. Therefore, 
this assumption fails. 

Let us consider another case, setting jRc = const at r ^ 1. Therefore, from the definition of 7(r) one can find 
that Rc{t) oc i^/^ and j{t) — > 0. Therefore, the term Ro^{u — 1) can be neglected at r ^ 1. Using the definition 
of the varia ble u one can set that the radius of inclusion increases in time according to the Lifshitz-Slyozov law 
R{t) cx 



46 



8 Simulations 



In order to verify our analytical results we performed numerical simulations by solving the system ( 19 1 in d = 2 



dimensional quadratic lattice of the linear size N = 256 with periodic boundary conditions with the mesh size 



£ = 0.5 and the time step At = 0.001. The noise term was modeled as the Ornshtein-Uhlenbeck process 47 
As initial conditions we have used (x(r,0)) = 0.5, (((5x)^) = 0.1. All measured quantities were averaged over 15 
independent runs. 

Typical patterns organized during the system evolution are shown in FigjS^ at a = 0.16. It is seen that at 
small times patterns of different sizes are formed. After the system selects patterns with a fixed size. At large time 
intervals spherical patterns are organized having preferably hexagonal symmetry as the Swift-Hohenberg operator 



predicts due to elastic interactions 23 24 30 . In Fig|8jD we plot the corresponding spherically averaged structure 
function at different times at relevant values for the wave-numbers to show the pattern selecting processes. It is seen 
that except the main peak related to period of main patterns there are well defined peaks at small wave-numbers 
(at elevated k the corresponding peaks have very small intensity, not shown here). During the system evolution its 
intensity decreases and only main peak is realized, its position is shifted toward 3?>f defining the period of stationary 
patterns. In insertion we have shown structure functions for two cases: a = 0.16 and a ~ 0.22. At a = 0.22 the 
period of patterns is smaller as the above stability analysis shows. 

Considering quasi-stationary limit (time intervals where {x), {{x — (x))'^) and {{x — (x))'^) are time-independent 
quantities, t ~ 10^) we can compute radius of inclusions by measuring surfaces of islands (number of points belonging 
to islands having spherical form). Evolution of the number of inclusion having fixed radius R at different a are 
shown in Figj9j It is seen that during the system evolution a distribution of islands becomes more compact and it 
has one the most probable value for R. 

After averaging this radius over independent runs we obtain the dependencies of the optimal characteristic 
length {Ropt){a) shown in FigjlOj where (• • • ) means average over experiments. On the top we present snapshots 
of patterns realized at different values for a at fixed e. It is well seen that graphs {Ropt){ce) repeat analytical 
dependencies Ro{a) in Figj?^ at small and large a where spherical islands forming hexagonal crystal structure are 
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Figure 9: Number of islands N versus the island size R at different times. Plots (a) and (b) relate to a = 0.16 and 
a = 0.3 Other parameters are: Dq = 3.0, rj = 0.3, A^'^^ = 0.1. 
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Figure 10: Averaged linear size of islands versus adsorption rate a at different e and Stationary patterns 

observed at e = 0.3 at different a are shown on the top. Other parameters are: Z?o — 3.0, tj — 0.3. 



possible. At intermediate values for a strip patterns are realized and no spherical islands are possible. Therefore, 
in Fig 10 there are gaps related to domains for a where no spherical islands are realized |^ Fitting procedure 
allows us to find that at small a where inclusions of the dense phase are possible the exponential law is realized, 
A — (Ropt) = B", where A and B are fitting parameters. At large a when islands of vacancies are realized one has 
exponential decaying, (Ropt) — A = exp{—a/B). Therefore, an introduction of the Swift-Hohenberg operator into 
the model for the system dynamics allows one to describe formation of lattice with different kind of inclusions: at 
small a one has molecular/atom inclusions forming hexagonal lattice, at large a one has a situation when deposition 
of atoms on metal surface can be described with formation of lattice of vacancy islands. 

To estimate the spatial size of patterns one can consider typical limits of rg ~ Inm and Ldif ~ Ifim. For the 
system of deposition of Al on TiN (100) the estimation in Ref. [14j shows that at room temperature with the lattice 
constant aAi — 4.05 x 10^^°m and pair interaction energy e = —0.22eV with lattice coordination number Z = A the 
size of patterns with ^ = Za\[e/T is around 30 x 10~^m. Therefore, patterns such as inclusions of dense/diluted 



phase are of nano-scales and are nano-patterns 15 . Patterns of diluted phase usually are realized as lattice of 
vacancy islands (for example in absorbed monoatomic layer of Ag deposited on Ru(OOOl) Ref. [?]). 

*In analytical study we have assumed that all inclusions are of spherical form 
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9 Conclusions 



We have studied pattern formation processes in stochastic reaction-diffusion systems with the relaxation flux obey- 
ing Maxwell-Cattaneo equation. The novelty of our approach is the usage of the phase field crystals formalism 
with special geometry of realized patterns and introduction of multiplicative noise obeying fluctuation-dissipation 
relation. 

Pattern formation is discussed on the model system describing adsorption-dcsorption processes. Considering 
pattern formation at early stages of the system evolution we have analyzed temporal stability of patterns and have 
shown that pattern selection processes are realized as a result of the diffusion flux relaxation to its stationary value. 

To study stationary patterns we have obtained the effective Fokkcr-Planck equation and its stationary solution. 
It was shown that multiplicative noise can induce non-equilibrium transitions at low energies for the adsorbate 
interactions. Using stability analysis we have computed diagram illustrating domains of the system parameters 
where stationary patterns can be realized. It was shown that introduccid multiplicative noise can induce re-entrant 
pattern- forming transitions where stable/unstable stationary patterns exist inside fixed intervals for lateral inter- 
action energy of the adsorbate. Stationary patterns organized during the system evolution form the crystalline 
structure of hexagonal symmetry. It is explained by introduction of the Swift-Hohenberg operator into the expan- 
sion for the interaction potential for the adsorbate. 

Considering quasi-stationary limit we have discussed evolution of radius of the dense phase inclusions into diluted 
one or vice versa. Using this alternative approach related to coalescence theory it was shown that all patterns are 
described by the characteristic radius Rq that less then the critical one Re- In the vicinity of Rc spatial perturbations 
grow exponentially, whereas in the vicinity of Rq such perturbations exponentially decay. We have found that both 
critical and characteristic radii have time asymptotics related to the Lifshits-Slyozov law. 
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